
pacman::p_load(tidyverse, stargazer)

load("Replication.Rdata")

# Table 3: Multirater inter-coder reliability ----
#Set up matrices for Krippendorf's alpha
alpha.matrix.del <- as.matrix(do.call(rbind, raw.icr.del))
alpha.matrix.disc <- as.matrix(do.call(rbind, raw.icr.disc))

#Krippendorff's alpha
alpha.del <- round(kripp.alpha(alpha.matrix.del, "ordinal")$value, 3)
alpha.disc <- round(kripp.alpha(alpha.matrix.disc, "ordinal")$value, 3)

#Set up matrices for Gwet's AC2 (missing values must be deleted)
gwet.matrix.del <- na.omit(as.matrix(do.call(cbind, raw.icr.del)))
gwet.matrix.disc <- na.omit(as.matrix(do.call(cbind, raw.icr.disc)))

#Gwet's AC2
gwet.del <- round(gwet.ac1.raw(gwet.matrix.del, "ordinal")$est[[4]], 3)
gwet.disc <- round(gwet.ac1.raw(gwet.matrix.disc, "ordinal")$est[[4]], 3)

t3 <- as.data.frame(rbind(c(alpha.del, alpha.disc), c(gwet.del, gwet.disc)))
colnames(t3) <- c("Delegation", "Discretion")
rownames(t3) <- c("Krippendorff's alpha", "Gwet's AC2")

stargazer(t3,
          summary = FALSE,
          header = FALSE,
          title = "Multi-Rater Inter-coder Reliability",
          type = 'text')

# Table 4: Use of flags by document type----
flag.tab <-
    data[, c("doc_type", "Ceremonial", "NotExecutive", "ConstraintsAppended", "Administrative", "Scored")] %>%
    group_by(doc_type) %>%
    mutate(n.ceremonial = sum(Ceremonial, na.rm = TRUE)) %>%
    mutate(n.notexec = sum(NotExecutive, na.rm = TRUE)) %>%
    mutate(n.const = sum(ConstraintsAppended, na.rm = TRUE)) %>%
    mutate(n.admin = sum(Administrative, na.rm = TRUE)) %>%
    mutate(n.scored = sum(Scored, na.rm = TRUE)) %>%
    select(c(doc_type, n.ceremonial, n.notexec, n.const, n.admin, n.scored)) %>%
    unique() %>%
    mutate(total.flag = n.ceremonial + n.notexec + n.const + n.admin) %>%
    mutate(total.docs = total.flag + n.scored) %>%
    t()


flag.tab <-
    flag.tab[c(2, 4, 5, 3, 6, 8), c(1, 2, 3)]

colnames(flag.tab) <- c("Executive Order", "Memo", "Proclamation")

row.names(flag.tab) <- c("Ceremonial", "Appended Constraints", "Administrative", "Not Executive", "Scored", "Total")

stargazer(flag.tab,
          summary = FALSE,
          header = FALSE,
          title = "Use of Flags by Document Type",
          type = 'text')

# Figure 1: Distribution of scores by document type ----

ggplot(data = data[data$Scored == 1,]) +
    geom_count(aes(x = Delegation, y = Discretion)) +
    theme_bw() +
    labs(size = "Frequency") +
    facet_wrap(~ doc_type)


# Figure 2: High-scoring documents over time ----
longdf <- gather(df, key = "measure", value = "value", c("yr.ct.big", "yr.prop.big"))

outcome.labs <- c("Number High-Scoring", "Proportion High-Scoring")
names(outcome.labs) <- c("yr.ct.big", "yr.prop.big")

ggplot(longdf, aes(x = year, y = value)) +
    geom_point() +
    geom_smooth(method = lm, color = "grey") +
    facet_wrap( ~ measure, ncol = 1, scales = "free", labeller = labeller(measure = outcome.labs)) +
    labs(x = "Year", y = "Number/Proportion of Documents") +
    theme_bw()



# Table 5: High-scoring documents over time ----

#Re-shape to year level
df_yr <-
    df %>%
    mutate(n = 1) %>%
    group_by(year) %>%
    mutate(n.yr = sum(n)) %>%
    select(c(year, yr.ct.big, yr.prop.big, n.yr))
df_yr <- df_yr[!duplicated(df_yr$year),]

#Re-scale document-level data year variable to start with 1
data_norm <-
    data %>% 
    mutate(year = year -1936,
           n = 1) %>%
    group_by(year) %>%
    mutate(n.yr = sum(n)) %>%
    ungroup()

#Baseline

m1 <- lm(yr.ct.big ~ year, data = df_yr)
m2 <- lm(yr.prop.big ~ year, data  = df_yr)
m3 <- lm(big ~ year, data = data_norm)
m4 <- lm(big ~ year*doc_type, data = data_norm)
m4b <- lm(big ~ year*doc_type + n.yr, data = data_norm)



stargazer(m1, m2, m3, m4, m4b,
          title = "High-scoring Documents Over Time",
          omit.stat = c("f", "ser", "rsq"),
          header = FALSE,
          label = "tab:over-time-models",
          dep.var.labels = c("Number High-Scoring", "Proportion High-Scoring", "High-Scoring"),
          model.names = FALSE,
          covariate.labels = c("Year", "Memo", "Proclamation", "Total", "Year x Memo", "Year x Proclamation", "Constant"),
          type = "text"
)

# Figure 3: Administrative delegation by document type ----
df_yr_type <-
    data %>%
    select(c(year, doc_type, big, hidel, hidis)) %>%
    mutate(n.total = 1)
df_yr_type <-
    aggregate(cbind(big, n.total, hidel, hidis) ~ year + doc_type, FUN = sum, data = df_yr_type) %>%
    rename(n.big = big) %>%
    mutate(p.big = n.big/n.total) %>%
    mutate(p.hidel = hidel/n.total) %>%
    mutate(p.hidis = hidis/n.total)

gg.grp <-
    pivot_longer(df_yr_type, c(n.big, p.big))

outcome.labs2 <- c("Number High-Scoring", "Proportion High-Scoring")
names(outcome.labs2) <- c("n.big", "p.big")

ggplot(gg.grp, aes(x = year, y = value, group = doc_type)) +
    geom_point(aes(shape = doc_type, color = doc_type), alpha = 0.5) +
    geom_smooth(method = lm, se = FALSE, aes(color = doc_type)) +
    facet_wrap( ~ name, ncol = 2, scales = "free", labeller = labeller(name = outcome.labs2)) +
    theme_bw() +
    scale_colour_grey() +
    labs(x = "Year", y = "Number/Proportion of Documents", color = "Document Type", shape = "Document Type")

# Figure 4: Administrative delegation by policy area ----
temp <-
    data %>%
    select(c(year, big, agg.topic, hidel, hidis)) %>%
    mutate(n = 1)

df.bytopic <-
    aggregate(cbind(big, n, hidel, hidis) ~ year + agg.topic, data = temp, FUN = sum) %>%
    rename(n.big = big,
           n.total = n) %>%
    mutate(p.big = n.big/n.total) %>%
    mutate(p.hidel = hidel/n.total) %>%
    mutate(p.hidis = hidis/n.total)

gg.grp.topic <-
    pivot_longer(df.bytopic, c(n.big, p.big))

ggplot(gg.grp.topic, aes(x = year, y = value, group = agg.topic)) +
    geom_point(aes(shape = agg.topic, color = agg.topic), alpha = 0.5) +
    geom_smooth(method = "lm", se = FALSE, aes(color = agg.topic)) +
    facet_wrap( ~ name, ncol = 2, scales = "free", labeller = labeller(name = outcome.labs2)) +
    theme_bw() +
    scale_colour_grey() +
    labs(x = "Year", y = "Number/Proportion of Documents", color = "Policy Area", shape = "Policy Area")

# Table 6: Trends by Topic
m1 <- lm(n.big ~ agg.topic, data = df.bytopic)
m2 <- lm(p.big ~ agg.topic, data = df.bytopic)
m3 <- lm(big ~ agg.topic, data = data)
m4 <- lm(big ~ agg.topic + year, data = data)

stargazer(m1, m2, m3, m4,
          title = "Trends by Topic",
          header = FALSE,
          dep.var.labels = c("Number High-Scoring", "Proportion High-Scoring", "High-Scoring"),
          covariate.labels = c("Foreign", "Year", "Constant"),
          omit.stat = c("ser", "rsq", "f"),
          type = 'text')
